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ABSTRACT 


Natural radio emissions from pulsating rad^o sources (pulsars) have 
been detected at Goldstone on a regular basis since 1968. Scientific 
analysis of these signals has stimulated ideas of traversing the inter- 
planetary medium and beyond, to the ’’home” of comets, using these natural 
beacons as navigation aids. Therefore, the techniques used in the ground- 
based observatiims of pulsars are described here in the required detail, 
many of them being directly applicable in a navigation scheme. The 
arrival times of *.he pulses intercepting Earth are measured at time 
intervals varying from a few days to a few months. Low-noise, wido-band 
receivers, unique to the stations of NASA's Deep Space Network, amplify 
signals intercepted by the 26-m, and 34-m, and 64-m antennas at Goldstone 
and in Spain. Digital recordings of total received signal power versus 
time are cross correlated with the appropriate pulse template, thereby 
providing an estimate of the pulse arrival* time relative to the station 
clock. Corrections are applied to the station clock to obtain arrival 
times relative to ephemeris time. Drifts in phase encountered during 
signal integration are removed. The arrival times are then referred to 
the bnrycenter of the solar system for scientific studies, and to the 
geocenter for export to other investigators. 
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SKCTION I 


INTRODUCTION 

Natural sources of pulsating radiation were first observed in 1967 
(Refs. 1, 2) by radio astronomers at Cambridge University in England. 
These pulsating radio sources (pulsars) are exciting because they (1) 
are thought to represent a stage in stellar evolution known as neutron 
stars (Ref. 3), providing a probe of an extremely dense state of matter; 
(2) exhibit changing, often dramatically so, behavior in the pulse per- 
iod, signifying changes in the spin rate of the neutron star driving 
the pulse generating mechanism; and (3) have large space velocities, 
implying violent origins. 

A series of measurements is being established using the NASA— 

JPL Deep Space Network (DSN) which will allow a probing of the neutron 
star interior and a direct measurement of changes in spin rate and 
angular motion. The purpose of the measurements is to measure the phase 
of the pulse train at known epochs. It requires little imagination to 
realize that the simultaneous comparison of the phase of a pulse train 
at two different DSN stations represents either a measurement of the 
difference between the station clocks, or a measurement of the loca**lon 
of one station relative to the other. In navigation, one assumes the 
clocks are synchronized. Differences in pulse phase are then attributed 
to differences in location. We need only to replace one DSN station 
with a spacecraft to complete the Illustration. 

The i>urpose of this report is to describe in appropriate detail 


Ll)c tfu'lmiqucs used currently in the DSN for measuring the arrival time 



of a particular pulse (al ternatively » the phase of the pulse train 
relative to a given epoch). Phase measurements were begun in September 
1968 and continue to the present. The 24 pulsars included in this series 
are listed in Table 1 with appropriates dates on Vvi.ich measurements began. 
An example of the results of this long series of measurements is presented 
in Table 7. 

Many characteristics of the pulsars have been or can be deduced 
from the measurement of pulse phase at a series of epochs if the epochs 
occur often enough (weekly to monthly"! and ovei a sufficient time span 
(years). Usually a parameterized model of the pulsar is constructed 
which is used to predict tne time-of-arr ival of a given pulse (or the 
phase of the pulse train at a given epoch). The parameters arc then 
adjusted to minimize the mean square difference between the predictions 
and the meas ^ments. The final values of the model parameters are 
dependent on the planetary ephemeris used in the data renu^'tion, since 
all observations are usually referred to the barycenter of the solar 
system. The importance of this is clear when one realizes that 2 milli- 
seconds flight time in the displacement of the barycenter amounts to as 
much as 1 arc-second displacement in the position of the pulsar. Errors 
in the knowledge of the barycenter with time scales of tens of years 
are caused by inaccuracies in the values of the masses of the outer 
planets. These slow curvilinear errors affect the values of the higher- 
order derivatives of the pulse period. Hence, though all current 
ephemerides may produce similar pulsai model paraimters, improveil 
ephetiieridcs civaiiahle twenty vears now will produce s ign i f i cant 1 y 

different results. Tables of geocentric arrival times, being important 
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Table 1. 

24 

of 

Pulsars Monitored at JPL 
2295 MHz and 2388 MHz 

at Frequencies 

Entry 


Pulsar 

Start Date 

1 


0031-07 

5 Mar 1970 

2 


0329+54 

5 Sep 1968 

3 


0355+54 

5 May 1973 

4 


0525+21 

1 Jun 1969 

5 


0628-28 

30 Dec 1969 

6 


0736-40 

1 Oct 1970 

7 


0823+26 

12 Feb 1969 

8 


0833-45 

22 Nov 1968 

9 


0950+08 

5 Sep 1968 

10 


1133+16 

25 Oct 1968 

11 


1237+25 

10 Aug 1969 

12 


1604-00 

1 Oct 1970 

13 


1642-03 

12 Jul 1969 

14 


1706-16 

1 Oct 1970 

15 


1749-28 

2 May 1969 

16 


1818-04 

29 Dec 1969 

17 


1911-04 

5 Mar 1970 

18 


1929+10 

29 Jun 1969 

19 


1933+16 

23 Dec 1968 

.M) 


2016428 

5 Sep 1968 

I 


2021+51 

30 Doc 1969 



2045-16 

10 Dec 1968 



21 1 1+^. 

30 Dec 1969 

J4 


22 1_7_+4 / , . _ 

1 ‘^69 


3 



archive's of pulsar bc'havior, are inciepondcnl of small f’han.;cs in l’im 
opheniGr ides of the planets. They form the exported fc^rm of tin; c^bser-* 
vational results. 

The measurements were performed earlier at a frequency of 2388 MHz 
on the 26-m antenna at Deep Space Station 13 (DSS 13) and the 64-m 
antenna at D5S 14. Current observations are at 2295 MHz, utilizing all 
facilities at Goldstone and in Spain. Detected pulsar signals are 
sampled, integrated and recorded on magnetic tape. Data collection and 
recording is discussed in Section II. The data are then edited and, for 
the early years, collated according to the pulsar, as described in Sec- 
tion III. The collated data, stored on magnetic tapes, was used to 
derive a template for each pulsar, representing the average power level 
across a pulse of radiation. The template is then cross correlated with 
each data record of the particular pulsar to obtain an estimate of the 
pulse arrival time. Tl.e procedures used in constructing the templates 
is described in Section IV and the cross correlation process is dis- 
cussed in Section V. The correction of the arrival time estimates for 
various small effects is discussed in Section VI. A sample list of 
geocentric tiraes-of-arr ival Is presented in Section VII. 


4 



4f 


I 

I- 

i 

£ 

i 

4 


I 


I 

I 





% 

f 

t 




.., , .if* it^-'^.. 


SECTION II 
DATA COLLECTION 

A. RECEIVING, DETECTION, AND RECORDING 

It was shown shortly after the discovery of pulsars that they could 
be detected at 2295 MHz (Ref. 4). Shortly thereafter, in August 1968, 

JPL observations at 2388 MHz were b' ;un using the 26-m antenna of DSS 13 
at the Goldstone complex, 'ilie wide bandwldths (about 12 MHz) and the 
low system temperatures (about l6 K) allowed detection of three of the 
first four pulsars discovered. Over the following months the observa- 
tions have been expanded to Include 24 objects, utilizing all antennas 
at Goldstone and Madrid. 

The receiver consists of the antenna, a right circularly polarized 
feed horn (left circularly polarized at 2388 MHz) , and a maser pre- 
amplifier followed by conversion to and amplification at intermediate 
frequencies (IF). Low-loss cabling then carries the signal to the 
receiver room for more IF amplification. I^e receiver bandwidth is 
limited to 12 MHz by either the IF portion of the receiver or by a 
filter Inserted for that purpose. Square-law detection of the signal 
is then performed, followed by low-pass filtering. The^ovtput of the 
low-pass filter is amplified by a wideband DC amplifier to a 1- to 
2-volt level. Subsequently, the signal is converted to digital form 
in a sampling process. The integration time (see Ref. 5 for a defini- 
tion) of the post-detection filter is chosen to match the time Interval 
between analog-to-dlgltal (A/D) conversions. ' 

Hie data processing, from sampling to the dumping onto magnetic 
tape, is under computer control. The detected signal is sampled at a 


4 


5 


4 





rate equal to 5000 times the apparent pulse period, corresponding to 
sample intervals of 50 psec up to HOO psec. Converted to a digital sig- 
nal by an A/D converter, the samples are stored in tho memory of an 
SDS 930 computer or a standard station >Iodcomp computer. Data samples 
are superimposed modulo the pulse period. Samples from 500 pulse periods 
are superimposed congruent ly. Then a 5000-word .array representing the 

superposition of 500 pulses is dumped onto tape. This accumulation of 

** - ■ . ^* % 

pulses continues until {a) detection is ensured, or (b) it is clear" the > 
pulses are too weak at this time to be observed. Accumulations rarely ., 
run more than one hour. 

The form of the recorded data is sketched in Figure 1. The pulse 
is sitting on top of a system noise power corresponding to a temperature 
T^. The peak equivalent temperature of the pulse is T^. The 5000 
samples spanning one pulse period are plotted from left to right. The 
time of the first data sample, plotted at the left, corresponds to a 
particular epoch of Universal Time (UT) as recorded at the station. 

This start time is recorded in UT (station). The train of received 
pulses then has a particular phase relative to the start time. The 
pulse train is folded in upon itself to beat down the effects of receiver 
noise fluctuations. Rather than think in terms of the phase of the 
pulse train, we ask for the arrival time of the first pulse after the 
start time. Hence, the pulse delay is added to the UT (station) epoch 
corresponding co the start time to yield an epoch for the arrival time 
of the pulse. 

The phase of the pulse train may not remain constant during the 
supr -position of many pulses. The subsequent drift in phase requires 
that a correction be applied to the measured pulse delay. This 




Figure 1. Nomenclature of the Pulsed Emission 



correction is discussed in Section Vf, D. Small corrections are required 
in the start time and in the UT (station) defined by a particular clock* 
These corrections are discussed in this section and in Section VI. 

B. THE .“:ONFTrjRATrON CODES 

The measurement system is shown diagrammatically in Figure 2* The 
experimental hardware is grouped into units defined by the dashed boxes. 
Each subsystem so defined has been interchanged with other similar sub- 
systems, causing significant changes in the measured phases. Hence, a 
system configuration code has been constructed to alert the analyst to 
changes in the system. The code varies from 0 to 21, representing 22 
significantly distinct combinations of antenna, receiver, observing 
frequency, computer, data collection controller, and station clock. A 
summary of the basic configurations (0-5) and (17-21) appears in 
Table 2. Additional configuration codes (6-16) were created to denote 
the use of an observing frequency of 2295 MHz or 2328 MHz, and their 
relationship to the basic codes are shown in Table 3. Left circular 
polarization (LCP) is to be assumed unless specified as right circular 
polarization (RCP) . 

Since the system configurations do not change in principle 
from one to the other, only configuration 0 is described in detail. 
Highlights of the other configurations are noted in the following para- 
graphs. Configuration 0 includes the 26-m antenna at DSS 13. The radio 
frequency bandpass is determined by the IF amplifiers. The final IF 
amplification, detection, and low-pass filtering is done by a receiver 
(labeled as receiver 1) built In 1968 by C. F. Foster. The timing con- 
trol box was designed by G. A. Morri*^ and is labeled timer I. 
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Figure 2. The Receiving, Detection, and Data Recording of Signals from Pulsars 




Table 2. Chronology of the Major Experiment Configurations for Collecting Pulsar Timing Data 






Table 3. The Configuration Codes Corresponding 
to Non-Standard Observing Frequencies 


Configuration 

Configuration 

Observing 
Frequency (MNz) 

Description 

M 

0 

2328 


9 } 


2295 


10 » 

1, 2, 3 

2328 


11 ) 


2295 


M 

4 

2328 

— 

il 


2295 


6 

5 

2295 





Denotes a correction w. s made to 

12 

Any 


standard post-detection Af^ 

13 

4 

2295 

RCP 

14 

5 

2295 

RCP 

15 

5 

2295 

DSS 13 via link to 




DSS 14 

16 

5i 

2295 

DSS 12 via link to 




DSS 14 

The relationship between 

the epoch-defining 1-second pulse and 

Universal Time 

as disseminated 

at the station, 

UT (station) , is shown 

in Figure 3 along with the relationships for the other codes. Exact 


relations between UT (station) and Coordinated Universal Time (UTC) Is 
discussed in Section VI, Data sampling in configuration codes less 
than 17 is performed by a 12-bit A/D converter with a 40 psec con- 
version time. The samples are stored in an SDS 930 computer having a 
16,384-word memory. The station telemetry processor assembly (TPA) 
is used in codes 17 or larger. 




(CHANGE IN WIDTH AS OF 1 NOV 1974) 









The final IF amplification and detection is performed in the steps 
shown in Figure 2. The first unit built was transported between DSS 13 
and DSS 14 until 10 Jan. 1974, when it was permanently installed at 
DSS 14. This unit was replaced. at DSS 13 by a programmable unit built 
by C. F. Foster (described in Ref. 6 and labeled receiver 2), This , 
change is not noted by a configuration code change. Hence, code 0 con- 
sists of an early phase (pre-10 Jan. 1974, Julian date 2442057.5) and 
a late phase (post-10 Jan. 1974). 

The timing control box used here is not described elsewhere, so 
a block diagram is presented here for tutorial purposes. Other control 
boxes used in other configurations operate on the same principle, but> 
differ in complexity and detail. Careful study of Figure 4 shows that 
a string of fast and slow pilses are generated by the timing control 
box to control the data sam])llng. The pulses are derived by counting 
down pulses from a frequency synthesizer, starting at a well-defined 
time by gating the synthesizer output. The gate must be armed and then 
a 1-sec pulse must be received before counting begins. To assist the 
operator in determining on which second the counting started, the l-sec 
pulses are counted down by 10 so that counting will start modulo 10 
seconds. 

The fast pulses, occurring every P /5000 seconds, where P is the 

i ® ® 

apparent pulse period, are used to strobe the A/D converter. The slower 
pulses, occurring every P seconds, are used to control the data transfer 
Into core memory. On the first slow pulse, a* computer interrupt is 
coincident with the slow pulse commands the A/D converter to begin 

$ 

j operation. About 20 psec later the end-of-convertion (EOC) Interrupt 

I forces the computer to accept the new word, and about 20 usee after 

I 

I 
( 


I 


13 










interruption, the word has been added to the proper memory location. 

Note that although it takes time to set up the computer to accept data, 
the EOC pulse arrives after the setup, so the first sample accepted is, 
in fact, coincident in time with the 1-sec pulse modulo one pulse period. 
However, a configuration exists in which the first sample is not coinci- 
dent, and a correction must be applied. 

The data sampling and superposition is under the control of a 
program written by G. A. Morris. The program also controlled the writ- 
ing of the results onto magnetic tape for further reduction. 

The apparent period P is the result of the Doppler effect caused 

a 

by Earth’s motion about the barycenter of the solar system. Tables of 

P are computed using a barycentric ephemeris of Earth and a crude but 
a 

adequate model of the pulsar. The model Includes the celestial 

coordinates of the object as well as the barycentric model of the pulse 

period (usually a period Pq at some epoch t^, and the first derivative 

Pn)* The values of P are tabulated at 15-mlnute Intervals between rise 
U a 

and set for any given day. The operator then chooses the value cor- 
responding to the middle of the intended Interval of pulse Integration. 

Configuration 1 consists of the 64-m antenna at DSS 14, the same 
IF amplification and detection unit, and the same timing box. used in 
the early configuration 0. A 12-MHz IF filter is added to the chain to 
limit the IF bandwidth. The computer Is again the SDS 930.- The A/D 
converter, however ^ responds within 5 psec of being strobed with a 
convert command. . computer , requiring about 20 psec to arm the 

IMP 

interrupt and start the data flow, misses the first data sample. The 
data stream is then shifted P./5000 seconds later in time. This shift 
is as large as 0.8 msec and Is added to the UT (station) recorded by 
the operator. 
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Timer 1 was designed to accept the rising edge of a positive 1-sec 
pulse as the definition of UT (station). The 1-sec pulse distributed 
at DSS 14 is negative-going, the rising edge being coincident with UT 
(station). The necessary pulse inversion then causes the timing box 
to begin counting earlier than intended by a time equal to the width of 
the 1-sec pulse. This amounts to 100 psec in configuration 1 and 
200 ysec in configurations 3 and 5. 

A programmable timing box, designed, built, and programmed by 
A. G. Slekys (Ref. 7) replaced the original timing box at DSS 14. This 
timer is labeled timer 2, and pertains to codes 2 and 3. A failure 
occurred in timer 2 which caused the start times of code 3 recorded 
between Julian dates 2442460.5 and 2443105.5 (17 Feb. 1975 and 
23 Nov. 1976) to be earlier than intended. The correction is a sub- 
traction of nearly 1 second. The exact correction was determined from 
a preliminary fit of a model of the pulsar emission times to the data. 

A correction of -0.996550 seconds to the start time aligned the affected 
arrival times with the times measured at DSS 13 (code 0 and 4) and with 
those unaffected from DSS 14. 

Timer 1 was replaced at DSS 13 on 1 April 1976 by a programmable 
data collector/timer designed by S. S. Brokl (Ref. 8) and programmed by 
K. I. ^k>yd. This timer is labeled timer 3, and pertains to configura- 
tion code 4. 

Occasionally a non-standard value is used for Af. in the post- 

Xt 

detection filtering. If a correction is applied, the configuration code 
is 12, corresponding to a correction to code 0 for dates prior to 1970 
and to code 5 after 1976, (None occur in the intervening years.) 


During the first half of 1979, the SDS 930 computer was phased out 
of service and removed in June 1979. This created th^ need to transmit 
IF signals over the microwave communication links to DSS 14 for data 
recording. The required 12 MHz bandwidth was obtained from the 6 MHz 
bandwidth of the link by folding the IF spectrum about the midpoint 
during base-banding. This technique was also used in obtaining signals 
from DSS 12. The configuration codes 15 and 16, then, indicate the 
recording at DSS 14 of signals from DSS 13 etid DSS 12, respidtiveryi '"*'' 
The most recent configuration for collectihg pulsar d^ta 
the use of the telemetry processor assembly (TPA) found in' All 
stations except DSS 13. The scheme was constructed by J; and M. Urech; ' 
at DSS 62. Many of the hardware functions of Figure 2 are duplioated 
by their software. The data output is still in the form of 5000H#ord 
records on magnetic tape. The configuration cod of this scheme for 
each station is listed in Table 2. 




EFFECTS OF RECEIVER PARAMETERS ON THE ARRIVAL TIME 


Pulsar signals are dispersive since the Interstellar plasma intro- 
duces a frequency-dependent group velocity. A broadband pulse, covering 
the radio spectrum from 30 MHz to 10 GHz, will arrive at the receiving " 
antenna and sweep through the receiver passband at a rate given by 

- KF^/DM (MHz/sec) (1> 

where the bandwidth B is in MHz, the observing frequency ~F is in MHz, 

-3 

At is in seconds, DM is the dispersion measure in parsecs (pc) cm , 
and K - 1.205043 x lO"^. Let Eq. (1) define the amount of smear At in 
the received pulse shape when all the other parameters are known. At 


F » 2388 MHz, a dispersion measure of 10 pr cm ^ causes thi% pulse to 





smear 75 ysec when the receiver passband Is 12-MHz wide. The smear is 
750 psec when the dispersion measure is 100. 

Computed values of smear are presented in Table 4. As long as the 
bandpass parameters remain fixed over the years, the dispersive effect 
should not affect the measured arrival times. However, changes in a 
filter element or the slight alteration in the tuning of the maser 
amplifier will change the frequency centroid of the passband. This, in 
turn, will change the time centroid of the detected pulse. The mag- 
nitude of this effect must be determined experimentally and has not 
been done. The efftict is expected to be small since it is not the 
limits of the passband that change, but only the shape. 

The response of the post-detection filter is chosen such that the 
integration time t is close to the time Interval P /5000, In the case 

3l 

of these simple RC filters, t = 2RC (Ref. 5), where R and C are the 
filter resistance and capacitance, respectively. Considerable care is 
exercised in using the appropriate filter combination for each pulsar 
throughout the observations. Use of non-standard values of t could cause 
up to 200 usee of time shift in the arrival time. The values are listed 
in Table 4. Note that the change from pre-10 Jan. 1974 (Code 0) to 
post-10 Jan. 1974 (Code 0) is principally a change in the post-detection 
filters, and much care was taken to duplicate the previous values of t 
in the new filters. 

A quantitative measure of the system performance is given by the 
ratio of o, the root-mean-square (rms) value of the noise fluctuations 
about the mean, divided by c« the mean. Letting N represent the 
number of pulses superimposed, 

T - ^ 

** /IJtB 
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Table 4. Pulse Smearing Due to Dispersion 
and Post-Detection Fil tering 


Pulsar 

Dispersion 
Ifeasure* 
(pc cm~3) 

Smear for 
B > 12 MHz 
(ysec) 

Integration 

Time 

(ysec) 

PSR 0031-07 

10.89 

79 

200 

0329+54 

26.776 

195 

50 

0355+54 

57.03 

416 

50 

0525+21 

50.955 

372 

400 

0628-28 

34.36 

251 

200 

0736-40 

161.0 

1174 

100 

0823+26 

19.4634 

142 

400 

0833-45 

69.08 

504 

50 

0950+08 

2.969 

22 

50 

1133+16 

4.8479 

35 

50 

1237+25 

9.296 

68 

200 

1604-00 

10.72 

78 

100 

1642-03 

35.71 

261 

100 

1706-16 

24.88 

182 

150 

1749-28 

50.88 

371 

400 

1818-04 

84.38 

616 

150 

1911-04 

89.43 

653 

150 

1929+10 

3.176 

23 

50 

1933+16 

158.53 

1157 

50 

2016+28 

14.176 

103 

50 

2021+51 

22.580 

165 

100 

2045-16 

11.51 

84 

200 

2111+46 

141.50 

1033 

200 

2217+47 

43.54 

318 

100 


*Values From Ref. 22 
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Oil'? expects for N = 500, t = 50 usee and B = 12 MHz, that o/\' ~ O.OU18. 
Throughout the years of data collection this parameter has varied between 
0.0013 and 0.0025, corresponding to 42 psec < t < 59 ysec. Hence, the 
values of Integration time listed in Table 2 were stable to within 
16 percent of the tabulated value. 
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SECTION III 


EDITING AND COLLATING 


A, EDITING THE DATA 

Successfully measuring the time-of -arrival of a pulse near a 
given epoch does not ensure that the result is free of biases. The 
editing process involves detecting the biases, correcting the measure- 
ment if the reason for the bias is known, or omitting the measurement 
if the reason for the bias is known but does not allow correction. A 
good example of a correctable situation is the use of a receiver fre- 
quency other than the 2388 MHz frequency used in most of the measurements. 

A preliminary model (position and pulse period) is devised for 
each pulsar. These models are used in a computer program of adequate 
accuracy to predict pulse arrival times at the barycenter of the solar 
system. A barycentric epheraeris of the Earth is used to refer the 
measured arrival times to the barycenter and a residual (measured minus 
predicted) arrival time is computed. The occurrence of a residual that 
is large compared to the average fluctuation is cause for further 
investigation. 

Several properties of the pulsar or the system are directly 
responsible for the need to correct or omit the arrival time measure- 
ment. They are: 

(1) Scintillations 

Interstellar clouds of ionized hydrogen impose variations 

on the reoelvod pulse power. The time scale of Che fluc- 

« 

tuations Is on the order of tens of minutes. If the pulsar 
Is far enough aw^^y, and several are, largo fluctuations in 
the pulse power force the signal down into the receiver 
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noise for many minutes, causing a spurious estimate of tlu- 
pulse arrival time. 

(2) Start Times 

The early timing box divides the incoming 1-sec pulse by 
10 (see Fig. 4), allowing the counting to be started at 
well-defined 10-sec intervals. It is imperative that the 
counter status register be synchronized with the station 
clock, modulo 10 seconds; i.e., on the even 10-second mark, 
the counter should read 0. Occasionally synchronization is 
lost. In these cases it is easy to deduce how many Integral 
seconds (never more than 10) the start time has to be 
retarded or advanced to obtain the true start time. 

(3) Receiver Frequency 

Dispersion caused by the interstellar plasma causes the 
arrival time to be a function of observing frequency. Using 
a non-standard receiver frequency produces a correctable 
error (see Section VI, F) . 

(4) Post-Detection Filter 

Occasionally a non-standard filter width was used, resulting 
in unusual time delay through the filter. A correction is 
readily applied. 

During the summer months of 1969 a bad detector connection 
caused unusually large (~1 ms) time delays. The corrections 
for this effect were derived by computing the autocorrelation 
function of the system noise, the width of which then gave 
the effective time constant of the filter. Both situations 
are labeled as system code 12 (see Section II, A). 
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(3) Miscellaneous 

In a very few cases the reasons for a large residual are 
unknovm. Common experimental situations In such cases 
are receiver tuning problems, repeated power failures, and 
radio frequency Interference (RFI). The measurement results 
are usually omitted. 

B. COLLATING THE DATA 

A data tape produced at Goldstone contains data files corre- 
sponding to the pulsars observed during a particular observing session. 
This format Is Inconvenient for producing average pulse profiles. (These 
profiles, discussed In Section IV, are a necessary part of optimally 
estimating the pulse arrival time). Therefore, the early data files of 
a particular pulsar were gathered onto one or more magnetic tapes In the 
order they were collected at Goldstone. Preliminary editing took place 
during this process, to delete obviously faulty data. The main cause 
for omission at this stage Is RFI. These data were used to compute 
the mean pulse profile presented In Section IV. 

Data collected after August of 1972 Is not collated according to 
pulsar. Data files are merely written onto magnetic tapes In the order 
they were collected without regard to the particular pulsar. The files 
were condensed by editing the data and filling each ta'pe completely. 
Estimates of arrival times are obtained after each tape Is filled (every 
2 to 3 months). 


23 




' -t 





SECTION IV 


PULSE-SHAPE TEMPLATES 

In the presence of white noise the maximum signal-to-noise (S/N) 
ratio for the output of a particular filter can be obtained when the 
filter response has the same form as the input signal (a matched filter). 
The required matched filter is realized by cross correlating a received 
pulse with an estimate of the noise-free shape of the pulse. The esti- 
mates of the noise-free pulse shapes are obtained by superimposing many 
thousands of pulses to average down the effects of noise. The result 
is referred to herein as a template. 

The proper registration or alignment of the pulses needs to be 
obtained to build a template with a minimum of smearing. Two passes, 
each described below, are used to optimize the registration. Correlation 
with a triangular function is first performed to obtain first-order 
estimates of pulse positicn within the data records. The template posi- 
tion in the data array corresponding to the maximum cross correlation 
coefficient is defined to be the pulse position. The numbers in each 
data record are then shifted circularly by the amount required for proper 
registration with all the other data records. Adding all the records 
produces a preliminary template. This template is cross-correlated with 
the original data to provide new estimates of pulse positions within 
the data arrays* The second result of the addition of all the properly 
registered data is taken as the correct average pulse shape. F .ally, 
the central 4096 words of the 5000 word array are used to form the tem- 
plate. See Section V regarding the choice of the 4096 words. 
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A. 


CORRELATION WITH A TRIANGLE 


In the early analyses of pulsar data (Ref. 9), the data were cross 
correlated with a triangular function closely matched in width to the 
pulsar signal. The pulses resemble the triangular function so the correla- 
tion function provides an adequate estimate of the pulse position within 
the data array. 

The shape of the triangular function lends itself to a particularly 
fast algorithm for computing the cross correlation coefficients. Let the 
data record be presented by f^, 1 ^ n £ 5000. Then correlate with the 
function 


Tr(n, n^, W) 



n-n 



W 


n-n^ 


< 1 


W 


(2) 


V 0 ; otherwise 

W is the half-width of the base of the triangle and n^ is the position 
of the peak in the data array. This function is shown in Fig. 5 
for Oq = W, W+1, and W+2. Inspecting Fig. 5, the correlation coefficients 
Cq and can be written as 


2N-1 

W-l 

2W-1 

■ 

■‘’X) w 


9 

II 

n»l 

n-Wfl 


2N 

w 

2W 


E 


~ ^ f + s ^ 

W n ^ 

2W+l-n , 
W n 

n»2 

n»2 

n-Wf2 



The difference between the coefficients becomes 
W 



n+1 


25 


( 3 ) 








POSITION 

Figure 5. Three Triangular Functions Overlying Data Array 


This is the first quantity computed. 

In a similar manner the difference between coefficients and 
ir 

W+1 

<=2 - ■ - 5 Z 

n=2 

This is the second quantity computed, but in the following manner. This 

new quantity is obtained from Eq. (3) by adding the term 

f 2f f 

_i W-1 ^2Wfl 

W " W " W 

In general the new difference in adjacent coefficients is computed from 
the previous difference by 

S ■ ^£-1 “ ^£-1 " ^£-2 


^ 5 »«-i 


^.^W+£-l ^2W+£-l^ 


(5) 
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The value of relative to the coefficient is calculated by accumulat- 
ing the sums ^i-l ~ ^£-2 through the many steps 

suggested by Eq. (5) : 







S- 


S-1 + ‘'i-l 


- C, = 


The position of the peak of the triangular function corresponding 
to the largest value of ~ denotes the position of the pulse in the 
data array. 

B. CORRELATION WITH THE TRUE TEMPLATE 


The cross correlation function C(t) of two continuous functions 
is written as 


c(t) = 



p(t)e(t - t) dt 


( 6 ) 


where p(t) is the pulse power as a function of time t, e(t) is the 
expected pulse shape (the pulse template) and i is the displacement in 
time of the two functions relative to some well-^def ined reference time. 
In the discrete case encountered here, the correlation coefficients can 
be evaluated by directly performing the operatic)!! suggested in Eq. (6). 
However, it is much faster to use Fourier analysis techniques. Direct 
evaluation of Eq. (6) requires 2,5 x 10^ multiplications, while using 
the fast Fourier transform (FFT) algorithm requires only about 1.4 x 10 


multiplications • 



Consider the Fourier tiansforms of the functions involved. The 


transform of the correlation function c(t) becomes 
C(s) = 

where s is the -transform variable. Interchanging the order of integration, 
C(s) = ^p(t) e du dx 

or 

C(s) - P(s) E(-s) (7) 

where P(s) and E(s) are the transforms of p(t) and e(t), respectively. 
Hence, the procedure is to calculate E(-s) for each template (calculated 
only once), P(s) for the data record, and then the complex function C(s). 
Inverse transformation of C(s) yields c(x). 

The FFT algorithm is used in a form that takes advantage of the fact 
that a Fourier transform of a real function is the sura of a real even 
function and an imaginary odd function. Hence, if the data array is L 
samples long, only L/2 values need be calculated in the transform. See 
Ref. 10 for the original work on this form of the algorithm or the detailed 
discussion of its use by Brigham (Ref. 11). Briefly, the data is placed 
in an array L samples long, where L is a multiple of 2. Pretending that 
the array is complex of dimension L/2, an FFT is performed on the complex 
array. A final step is needed to unscramble the result into (L/2) - 1 
complex pairs, one 0-frequency value, and one Nyquist-rate value. Alter- 
natively, putting complex data such as the quantities P(s)E(-s) in this 
form, the Inverse of the above process can be formed. Software to per- 
form the calculations was written for this project by G. A. Morris. 


I 


p(t) e(t - x) 


j.. -i2irsx . 
dt e dx 
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The central 4096 data points of the integrated pulse profiles were 

selected for the templates since the correlation process Involves using 

12 

the FFT algorithm for 2 data points. The correlation 1s done twice 
to cover all 3000 data points. The details are presented In Section V. 
The value of t corresponding to the maximum value of c(t) denotes the 
position of the pulse In the data array. 

C. SUBTRACTION OF THE BACKGROUND 


Once an estimate of the pulse phase has been made, the pulsed rjwer < 
must be separated from the steady component contributed by the system. 

Let the signal be represented by a pulse template U(t - t) of unity maxi- 
mum height. The pulse is superimposed on a constant term N due to system 

s 

noise. Associated with are small noise variations n(t) caused by a 
non-zero post-detection bandwidth. Hence, the total signal S(t) becomes 


S(t) = gU(t - t) + n(t) + N (8) 

s 

where g scales the template to the data. The expected value of n, 
denoted by E{n}, is zero. The objective is to minimize the quantity 


■/ [ 


S(t) - gU(t - t) 




dt 


(9) 


All terms in are constant except the term 


2g 



S(t) U(t - t) dt 


which must be maximized to minimize But this is merely the cross 

correlation performed above. 

Estimates of g and N are needed to separate the pulse from 

s 

receiver noise. The estimates should minimize 



I 


S(t) U(t - t) dt + g / U (t - T) dt + 


/ 


-./• 


(t-T) dt 
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In 


i . 0 • - J St. t) dc + g fvtt - t) dt + 


»s/ 


dt 


where t is the estimate of x from the correlation procedures. Solving 

for the estimates g and H of g and N , respectively. 

s s 


/"/ 


- J s(t) dt J 


S(t) U(t - t) dt - / S(t) dt / U(t - t) dt 


/D 


N 


Jsit)dt j U^(t-x) dt- J 


(t-x) dt- / S(t)U(T-x)dt /U(t-x) dt 


/■ 


/D 


( 10 ) 


where 


D = /dt/«= 


(t - x) dt - 


/ 


U(t - t) dt 


Note that only the quantities involving S(t) need to be calculated for 

each data record. All other terms are calculated once for each template. 

The discrete version of Eq. (10) is used in the data reduction. 

The value N is subtracted from each data record before adding the 
s 

record to the evolving template. 

D. THE TEMPLATES 

Several hundred thousand pulses were superimposed as described 
above to form the final estimates of the average pulse shapes of 24 pul- 
sars. The results appear In Figures 6 through 29. Part (a) of each 
figure represents the template used in the data reduction. The abscissa 
presents time in terms of the pulse period P, while the ordinate Is In 
arbitrary units of power. The vertical marks above and below the main 
pulse define the time origin. Some templates were smoothed beyond the 
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amount inJierent in the receiver. The smoothing was obtained by convolving 
a sinc(c) function of the appropriate width with the original data. The 
amount uf smoothing noted in the figure caption represents the interval 
from the peak to the first zero of sinc(t). A quantitative discussion 
of the choice of smoothing is presented in the next section. 

Part (b) of each figure presents the details of the pulse on an 
expanded time scale. If the data has been smoothed more than in part (a) 
of the figure, the amount of smoothing is noted. The position of the 
pulse relative to the vertical marks above and below the main pulse 
defines, for that pulsar, the "zero phase" condition relative to a given 
epoch; l.e., to obtain a zero time-of-arrival residual relative to a 
given epoch, the vertical marks must be coincident with that epoch. 

E. TEMPLATE SMOOTHING 

The accuracy with which one can estimate the arrival time of a pulse 
Increases as the random fluctuations in the signal decrease and as the 
sharpness of the pulse inc^'eases. However, attempts to smooth the signal 
to lower the effects of receiver noise also tend to destroy the sharpness 
of the pulse. A compromise must then be reached. This is done by esti- 
mating the uncertainty in the arrival time measurement when the signal-to- 
nolse ratio is large. 

Imagine that having performed the correlation process to measure 

A 

the arrival time of the pulse, an estimate t has been obtained. The 

A 

scaling of the template, g, and the constant component, N , have also 
been estimated. Allowing for small errors in the estimates, the model 
of Eq. (8) Can be rewritten as a Taylor expansion: 

S(t) » gU(Ax) + AgU(Ax) - gAxU'(Ax) + N„ + AN„ + n(Ax) (11) 

S 5 






I 

Figure 6. The template of PSR 0031-07: (a) P = 0.943 sec; 

(b) Expanded time scale 



Figure 7, The template of PSR 0329+54: (a) P ■ 0.715 secj 

(b) Expanded time scale 
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Figure 8. The template of PSR 0355+54; (a) P « 0.156 sec. The pulse 

has been smoothed by a window 160 psec wide, (b) Expanded 
time scale 



Figure 9. The template of PSR 0525+21. (a) P « 3.745 sec. The pulse 

has been smoothed by a window 2.2 ms wide, (b) Expanded 
time scale 
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Fxgure 10. The template of PSR 0628-28; 
(b) Expanded time scale 


(a) P = 1.244 sec; 



Figure 11. The template of PSR 0736-40: (a) P 

(b) Expanded time scale 


0.374 sec; 
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Figure 16. The tenplate of PSR 1237+25: (a) P = 1.382 sec; 

(b) Expanded time scale 


% 



Figure 17. The template of PSR 160A-00: (a) P = 0.422 sec; 

(b) Expanded time scale 
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Figure 18. The template of PSR 1647.-03: (a) P = 0.388 sec; 

The data has been smoothed by a window 230 psec 
wide, (b) Expanded time scale 



Figure 19. The template of PSR 1706-16; (a) P = 0.653 sec; 

The data has been smoothed by a window 390 psec 
wide. (b) Expanded time scale 










(a) 



20 ms 


Figure 21. The template of PSR 1818-04; 
(b) Expanded time scale 


(a) P = 0.598 sec; 
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Figure 24. The template of PSR 1933+16: (a) P = 0.359 sec; 

(b) Expanded time scale 



Figure 25. The template of PSR 2016+28; (a) P - 0.559 sec; 

(b) Expanded time scale 
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Figure 26. The template of PSR 2021+51: (a) P ® 0.529 sec; 

(b) Expanded time scale 
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Figure 27. The template of PSR 2045-16; (a) P « 1.96^ sec. 

The data has been smoothed by a window 1.2 ms wide, 
(b) Expanded time scale 


t 
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Figure 28. The temqlate of PSR 2111+46; (a) P = 1.015 sec; 

(b) Expanded time scale, with the data smoothed 
for display by a window 2.0 ms wide 


I 



Figure 29, The template of PSR 2217+47: (a) P = 0,538 sec; 

(b) Expanded time scale, with the data smoothed 
for display by a window 1 nu> wide 
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where the prime denotes a derivative, and 


Ag = g - g, AN = Ng - Ns and At = t - T 

s ^ 

Substituting Eq. (11) into Eq. (9) yields the mean square error 

^ ^ /N 

in using estimates g, N and T instead of the mean values g, N and t. 

s s 

Minimize the mean square error with respect to At to obtain 




AgU(Ax) - gAtUMAi) + AN_ + 


n(Aa)J 


[AgU‘ 


(At) - gU*(Ax) - gATU**(AT) + n*(Ax) I dx = 0 


n' (A t)J 


where 


U’ * dU/dt, U" = d^U/dt^, etc. The integration is over one full 
pulse period. Collecting terms, we make use of the following facts: 

(1) U(0) = U(P) * 0 

(2) U'(0) = U’(P) = 0 
dx = n(0) - n(P) « 0 


and 


Then 


(3) 


(4) 




dx 0 


h 


dt 


Ax 


/ 


g / (U’)^ dx 


where terms of the order Ag/g have been omitted. 
Since 


E{Ax} 


fv EMir/i / dx . 3. 


( 12 ) 



then 


0 ^ " E{At^} = 


I jvHs) nW ds / U'(q) n 

J (U')^ dT 


(q)dq 


u’(s) U'(q) Rjj(s- q) ds dq 


J (U')2 dx 


where R^(At) Is the auto-correlatlon function of the random fluctuation 
2 

n{t). R (0) = a . If R always decreases significantly before U* can 
n n n 

change significantly, then this result remains independent of the post- 
detection bandwidth. Measuring the variance of the fluctuations In a 
region not containing the pulse, we can write 


’'“f- Tt 

j (U')^ dx 
. 0 

J 


The denominator is corrected for contributions to U' from random fluc- 
tuations. The expected standard deviation a? jclated with x is then o . 

X 

The expression for satisfies our intuition with respect to how 
well X can be measured for a given slgnal-to-noise ratio (g/Ojj) and a 


given pulse sharpness (measured by 


/ 


) dx). A test of Eq. (13) was 


performed which dio not rely on intuition. The* extremely noise-free 
template for PSR 0833-45 was used as the basis for producing 100 records 
of simulated data. A random number generator was used to superimpose 
uncorrelated fluctuations onto the original template. Equation (13) was 
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used to choose the variance a given that o = 5 Msec. The sharpness 

f 2 

I (U') di is known from the template, and g = 1 in this case. 

Would the correlation method define a random process t from the 
100 simulated data records wlch the correct standard deviation? It did. 
The true value of t was known in this test, and the standard deviation 
o£ the estimates t about the expected value t was 5.2 psec. 

The templates were smoothed to the resolution where o^ stopped 
decreasing due to a lower and began increasing due to loss in sharpness 
of the pulse. 

F. CENTRAL SPIKES 

Occasionally, the pulse strength in some of the weaker pulsars 
ebbs to the point of invisibility (a S/N ratio of 1 or less), even aft®r 
an hour of integration. In these data records the template seeks out 
the stronger, sharper random fluctuation... The resulting estimate of t 
represents the position in the data array of one of the larger noise 
spikes. In superimposing data records to produce a template, a narrow 
central spike was formed. 

The central spike phenomenon was observed in the templates for 
PSR 1818-04 and PSR 2111+46. The spike was removed by placing in the 
central nine array locations a random process with a mean computed from 
the eight neighboring values. The variance was set equal to that of the 
pulse-free portion of the array. 

G. DEFINING THE ZERO-LEVEL 

Defining the zero-level of the pulsed radiation is straightforward 
in the cases where the pul . js are narrow and confined to one window. 

This is not to say that the radiation between the pulses is not in part 
pulsed and associated with the pulsar. In most cases, however, the 
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emission is featureless. If a portion of the flat interpulse region 
is also pulsed, measurements of a different nature are required to detect 
this possible situation. In general, the mean value of the interrulse 
emission is used to define the zero-level. 

In several cases, such as PSR 095&f08, the presence of sharp or 
broad features requires a careful definition of the interpulse region. 

In these carefully defined regions the mean value is then taken to be 
the zero-level of pulsed emission. 



SECTION V 


ESTIMATING THE PULSE ARRIVAL TIME 

Each data record on magnetic tape presents an opportunity to 

measure the arrival time of a received pulse relative to a given epoch. 

The measurement is completed by correlating the appropriate template 

(Figures 6 through 29) with data records represented by Figure 1. The 

correlation process leads to an estimate of the maximum cross-correlation 

coefficient. However, the coefficients occur only every P /5000 seconds 

a 

of delay, so interpolation in the neighborhood of the largest coefficient 
is required to obtain the best estimate of the delay corresponding to the 
true maximum of the coefficients. The value of the largest correlation 
coefficient and an estimate of the random fluctuations due to receiver 
noise are ured to estimate the uncertainty in the time delay. 

A. THE CORRELATION PROCESS 

The data records are 5000 words long. Two arrays of 4096 words 
each are ccnstrucced to allow easy use of the FFT algorithm. The two 
separate results, each 4096 words, are combined to yield one 5000-word 
array of cross-correlation coefficients. 

The template is shifted in position before transformation, as 
shown in Figure 30. The time within each pulse corresponding to the 
ar rival time of the pulse is, in practice, the time marked **0** in Fig- 
ures 6 through 29. Each template is shifted circularly until this time 
occurs in the first word of each array. Hence, the position of the 
largest correlation coefficient vlll be the estimate of the arrival time 
of that particular pulse. In addUion, this shift keeps the size of 
the imaginary component of the transformed template to a minimum. 


TEMPLATE 

ARRAY 



TRANSFORMED 

ARRAY 



REAL 


IMAGINARY 


Figure 30. Shifting and Transforming the Template 
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obtained when no smearing is present represents a bias that can be 
estimated and removed. The bias is estimated by first calculating the 
expected arrival time based on the model of the pulsar, and then compar- 
ing that with the arrival time expected for a pulsar with a constant 

/ 

period and constant Doppler shift. The mean of the differences is the 
estimate of the bias. 

The expected arrival time t^^ of the pulse arriving after the 
starting epoch is given by 

‘k- "“'o'] 

where is the elapsed time between k pulses at the barycenter of the 
solar system, d(t) is the distance of the received wavefront from the 
barycenter, and c is the velocity of light. Expanding the barycentric 
pulse period P(t) in a Taylor series about t„, 

\J 

It can then be shown that 

The quantity d(t) can be replaced with the scalar product -r (t) 

a 

•Up(t), where r^(t) is the position vector of the antenna relative to 

the barycenter at time t. The unit vector U is the direction vector 

P 

of the pulsar as seen from the barycenter. Referring to Figure 34, note 

that |r I » a 1 AU, and R/a >> 1. Using the facts that 
a 

Y ■ -a cos 0 - d 

2(R-d)Y - 
and 

+ Y^ - Z - d^ + a“ + 2ad cos 6, 
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Figure 31. Computing the Cross Correlation Coefficients 
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Figure 32. An Example of the Cross Correlation Process; 

PSR 1818-04 on 27 Jan. 1970 at 16'> 58'" 00^ UT 

B. ESTIMATING THE DELAY 

The correlation coefficients are examined to determine the position 
in the array of the largest coefficient. In the cases of PSR 0355+54 and 
PSR 0833-45, two and three pulses, respectively, occur per data record. 

In these cases two or three positions are determined, eacli separated by 
a pulse period, and averaged together at the end of the estimation 
process. 

Subsequent refinement of the estimate of the pulse position is 
obtained through interpolation. Interpolation is obtained by fitting, 
in a least-square sense, a fourth-order polynomial to M coef f icitMits 
either side of the largest coefficient. At position (measured rela- 
tive to the position of the largest coefficient) in the array of' 
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correlation coefficients, a predicted value of a^ is obtained. 


1“0 


Note that x , - x = 1. We have for the 2M + 1 coefficients Involved, 
i+X i 


^ M “ u + ••• 3/X « 

-M 0 1 -M 4 -M 


^0 " ^0 


C- ** ••• 3/X- 

1 0 11 4 1 


St “ + ^l^M + • • • ®aSi 


The polynomial coefficients Bj need to be determined* Denoting the meas- 

ured coefficient by C^, we want to minimize the quantity 

2 


M 


' L ''i - ■ E 


i=-M 


- E 


1=0 


3S 

Requiring -r — = 0 for i * 0 to 4, we find the set of equations 


M 






i=-M 

M 


= 1; “ ^0 2 ^ ®1 Z] "'l ®4 Z^ ’^i 


(14) 


1=-M 


1 = 4 


= E ^1*1 ■ °o S *1 ^ ‘i 2 *i * ••• *\ Xj *i 


i=-M 
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Since the occur at uniform intervalr, the sums of odd powers of are 
zero, and the sums of even powers of x^^ are precomputed for each pulsar. 
The size of M is chosen separately for each pulsar. Within the range 
x_j^ < < Xy, the fourth-order polynomial provides a good fit to the 

correlation coefficients obtained when the template is correlated with 
itself. 

The system of equations in Eq. (14) is solved for the polynomial 

coefficients a^. The curvature of the resulting polynomial is computed. 

A positive second derivative indicates a distortion due to a low signal- 

to-noise ratio, and the Interpolation is not attempted. Otherwise, the 

position Xp at which the polynomial peaked is determined. Then, stepping 

in intervals of 0.01 between x -1 and x +1, the refined estimate x. of 

P P d 

where the polynomial reaches its peak is determined. 

The arrival time t, of the first pulse occuri.ig after the start 

time t is readily calculated. The separation between data array eleme*.cs 
s 

is P /5000 seconds. Hence, t = t + x , P /5000. 
a ’ a s d a 

The interpolation performed near the peak coefficient in Figure 32 
is presented in Figure 33. The measured values of the cross-correlation 
coefficients (solid circles) are displayed versus the corresponding array 
location, which is proportional to time delay. The solid curve is the 
polynomial resulting from the least-squares fit. The pea, of the poly- 
nomial is near x » 589.0, The region near this value of x is enclosed 
by the broken box, and reproduced on a finer scale in the inset. Pro- 
gressing from X = 589 to x * 590 in steps of 0.01, the peak, denoted 
by the vertical arrow, is found at * 589.79. 
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Figure 33. Correlation Coefficients Versus Time Delay 
Near the Peak Shown in Figure 32 


C. ADJUSTMENTS 

If the template has been shifted due to a large asymetry about 
the zero-phase position, the effect of that shift is removed at this 
point in the data reduction. If two or three estimates of were made, 

they are averaged together at this point. 

D. UNCERTAINTY IN THE ARRIVAL TIME 

Each estimate of the arrival time t ii: accompanied by an estimate 
of the uncertainty o^. A modification of Eq. (13) is used. The scale 
factor g and the baseline level N are first estimated using Eq. (10), 
Quantities dependent on the template U(t) only, having been computed 
earlier In the template generating process, are on magnetic tape i’Uh 
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the template itself. Note that the quantity S(t)U(t - t) dt is merely 
the maximum correlation coefficient. The root-mean-square (rms) noise 
level includes a contribution by the template. The denominator of 
Eq. (13) is another quantity computed earlier. 
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SECTION VI 

GEOCENTRIC ARRIVAL TIMES 

Now that the arrival time of the pulse has been measured relative to 
an epoch defined by a particular clock. It Is time to consider referring 
these station arrival times to the geocenter. There are several correc- 
tions which must be applied to the measured arrival time before useful 
geocentric results are obtained: (a) A bias In the arrival time, caused 

by the smearing of the pulse during superposition, must be estimated and 
removed; (b) the geometric correction, changing the spatial reference 
from the antenna to the geocenter, must be applied; (c) an extra delay 
J must be removed which Is caused by using frequencies lower than 2388 MHz; 

(d) a small annual effect caused by the Interaction of the Doppler and 
dispersion phenomena must be removed; (e) the clock epochs used In the 
measurements must be expressed as Coordinated Universal Time (U'lC) epochs; 
m (f) UTC epochs must then be expressed as ephemeris time (ET) epochs. The 

result is then the arrival time that would have been obtained had the 
measurement been performed at the geocenter using a clock running on 
ephemeris time. 

A. CORRECTION FOR PULSE SMEARING 

A small amount of smearing Is inevitable as the received pulses 
are superimposed. The Doppler shift caused by Earth's rotation and 
orbital motion changes over the interval of integration. Also, differ- 
ences can exist between the actual pulse period and the assumed period. 

The result is a slow drift in the pulse phase -elative to the starting 
epoch, causing the centroid of the Integrated signal to drift forward 
or backward In time as pulses are superimposed. The difference between 
the pliasc of the centroid of the superimposed pulses and the phase 
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obtained when no smearing is present represents a bias that can be 
estimated and removed. The bias is estimated by first calculating the 
expected arrival time based on the model of the pulsar, and then compar- 
ing that with the arrival time expected for a pulsar with a constant 
period and constant Doppler shift. The mean of the differences is the 
estimate of the bias. 

The expected arrival time t^^ of the pulse arriving after the 
starting epoch t^ is given by 

where t^^^ is the elapsed time between k pulses at the barycenter of the 
solar system, d(t) is the distance of the received wavefront from the 
barycenter, and c is the velocity of light. Expanding the barycentric 
pulse period P(t) in a Taylor series about tg, 

p = Pq + ^0^*^ ■ V 

It can then be shown that 

The quantity d(t) can be replaced with the scalar product -r (t) 

a 

•Up(t), where r^(t) is the position vector of the antenna relative to 

the barycenter at time t. The unit vector U is the direction vector 

P 

of the pulsar as seen from the barycenter. Referring to Figure 34, note 

that |r I « a 1 AU, and R/a >> 1. Using the facts that 
a 

Y * -a cos 0 - d 

2(R-d)Y - 
and 

2 2 2 2 '* 

X + Y - Z - d + a“ + 2ad cos 0, 
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Figure 34. Geometry for Calculating the Phase Drift 


the distance d becomes ^ 

d s -a cos ® ~ 

where R/a » 1. (The smallest value of this ratio is about 10^.) 

Consider the arrival time obtained if P(t) - P^ over the entire 
interval of Integration. Suppose there is a delay of I pulses from the 
time the timing pulses are started to the time when the data collection 
is started (usually a few seconds at most) . Then M groups of N pulses 
(usually 500) are superimposed, with n pulses elapsing between groups 
to allow for a visual display of the progress of the observation. Fig- 
ure 35 displays the numerology used in counting these pulses. Under the 


! 
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Figure 35. Nomenclature for Enumerating Pulses 

assumption of a constant period, the arrival time t^j^ of the pulse 
after the start epoch tQ Is 


taj^ = + j + (N + n)(m - 1)] = kP^ 

th 

where j is the pulse in question in the m block of N pulses, and P 

a 

is the assumed constant period. (The synthesizer frequency f used in 

s 

the observation, as per Figure 2, Is 10^/P .) 

a 

The mean of the difference tj^ - t^j^ Is given by 

M N 

® {^k ■ ’^ak} “ ‘^k " *^ak “ Xrf £^*^k " 

m-1 j-1 


Let D • tj^ - t^j^, representing the bias caused by the smearing. Sub- 
stituting Eqs. (15) and (16) Into the expression for D yields 


^^0 " 




M 


NM 


m*l 


N 

J'l 


p.p 


+ «, + (N + n)(m - 1)} + 


M N 

0 0 ^ 

imZrnJ /mJ 

m=l j=i 


(17) 


M N 
m"l J«1 


(d(tj^) - d(tQ)} 
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The mean of the distances can be replaced by the distance near the 

midpoint of the range of j, yielding for the third term of Eq^ (17), 


M 


- V 

cM 


m=l 


N+1 


+ i + (N + n)(m - 1) 


- 'a<V 


• M (t ) 


Performing the indicated summations, Eq. (17) becomes 

(Pq-P^) 

+ !A j + 4 + (N nHM ^ (N^ - 1) H- n)^ -...U j 

M 

z 


(18) 


+ 4 + (N + n) (m - 1) 


- 


%<‘0> 


It has been assumed that the signal strength is constant over the range 
of k. Large fluctuations cause Eq. (18) to be in error, and if the drift 
is large, the results may be rejected. 

The expression derived above for the drift is applied tc all the 
data. In all cases, N > 500. The delay I is usually 0. The plot time 
n is equivalent to 30 seconds before Julian date 2440195.5, and 17 sec- 
onds thereafter. The numbe ' M of 500-pulse integrations varies depending 
on pulse strength, ranging from the equivalent of a few minutes to over 
two hours. 

The Jet Propulsion Laboratory ephemerlr OE 96 is used to compute 
bary entric positions r . One major export from this observing program 
is ephemetls Independent, geocentric arrival times. However, the drift 
correction is needed, and is best done by the observer early in the 
reduction process. The use of a different ephemeris of similar quality 
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will not change the corrections of Eq. (18) significantly. A large drift 
correction is on the order of a few milliseconds, and a velocity differ- 
ence between two ephemerides of greater than 0.1 percent would be required 
to produce a 1 usee change in the drift estimate. Ais is usually less 
than the error estimates encountered in these measurements. Note also 
that this correction is dependent on a good knowledge of the pulse period 
and its derivatives. These have been well-determined in a separate data 
analysis effort. 

The sum of the barycentric position vector of the geocenter and 

the geocentric position of the antenna, r^, yields r^. The computation 

of r is discussed in the next section. The direction vector U is 
g P 

derived from a model-fitting procedure, which is beyond the scope of 
this report. All the pulsar positions used here are known to at least 
1 arc-sec, an error that will produce a geocentric arrival time error 
of 0.1 usee or less. 

An example of phase drift and the subsequent correction for the 
effect is shown in Figure 36. Sixteen groups of 500 pulses of PSR 0833**45 
were examined separately to estimate x. The solid circles represent these 
estimates relative to the estimate for m - 1. The line through the solid 
circles represents a fit-by-eye to the linear drift. The triangle is 
the result of examining the accumulation of the 16 blocks of data. As 
one would expect in a well-behaved situation, this estimate of x has 
drifted about one-half the distan<^ ^rom ”0” as the estimate for the 
m = 16 case alone. This line through the triangle and the zero for 
this example (0 at m = 1) delineates the hypothetical drift one would 
find in examining ar ’ruulation of the preceding blocks of data. The 
quantity D, derived from Cq, (18), is the amount the result represented 
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Figure 36. An Kxample of Drift in the Phase of PSR 0833-A5 
on 1 Dec. 1968 at 00 ^ 00« (UT) 

by the triangle is siiifted to correct for drift. Note that by shifting 
the datum down by D ms» it is nearly in line with the intersection of 
the dashed l?ne and the m = 0 axis. In a later analysis the resulting t 
was found ae only 15 psec away from the arrival time predicted by a 
smooth model, which was fitted to several months of data. 

B. ANTENNA POSITION 


The referral of the measured arrival time to the geocenter requires 

the addition of the time delay r • U /c, where r is the geocentric 

Pg g 

position of the antenna and U is the geocentric direction vector of 

P8 

the pulsar. In practice U The coordinate system used here is 

Ph P 

the right-handed rectangular system (X, Z). The coordinate X is in 

the direction of 0^^ Kipju Ascension (FA) at Tlu' station position 

vit a pa r l i lU I .• I T I'T t fun'i) is calculated from llu oujto ua coi»rd i aa t o s i o. 




Table 5 and the knowledge of the Greenwich Hour Angle (GHA) • The station 
coordinates are a version of Location Set (LS) 44, which are based on the 


cphemeris DE 96 (Ref. 19). Tlie GHA is computed as outlined in Table 
(see Ref. 20). 

Table 5. Geocentric Antenna Coordinates 


Antenna 

West 

Longitude 3 
(deg) 

Latitude A 
(deg) 

Radius R 
(km) 

DSS 11 (26m) 

116.849390 

35,208047 

6372.011 

DSS 12 (34m) 

116.805462 

35.118665 

6371.994 

DSS 13 (26m) 

116.794863 

35.066546 

6372.117 

DSS 14 (64m) 

116.889507 

35.170847 

6366.227 

DSS 62 (26m) 

4.367811 

40.263145 

6369,963 

DSS 63 (64m) 

4.247991 

40.241318 

6370.048 


Table 6. Relations Used in Computing the GHA of an 
Observing Station 


At the epoch d + SEC/86400, 

c:HA = 100.0755426 + d(0. 985647346 + 2.9015 x 10~^^ d) 
+ 0 SEC + fiy cos(u) degrees 


where: (1) 

( 2 ) 

(3) 

(4) 

(5) 


d = Julian days pasl 2433282.5 
SEC = seconds past 0 on day d 

6 = 4.17807417 X 10~^/(1 + 5.21 x lO"^^ d) deg/sec 

6Y, the nutation in longitude, is taken from 
ephcmeris DE 96 

e = 4.09205684 x 10~^ 

+ T(-2. 28534 x lO"^ + T(-l. 5 x 10“® + 8. 7 x lO"^)) 

where T Is in Julian centuries of 36525 days past 
2433282.5. 


64 


radians 










The components of the geocentric antenna vector are computed in the 
1950.0 coordinate system, requiring rotations for precession and nutation. 


Considering an antenna at longitude 0 and latitude A at a given epoch. 


r = 
g 


r 

gx 


■"ii 

"12 

**13 


1 

-6ycos(e) 

-6ysin(e) 

r 

gy 

= 

^21 

P22 

**23 

X 

Aycos (T) 

1 

-6e 

r 

L 8z. 


/ 3 I 

P32 

^33 


Ay sin (7) 

6e 

1 


gx 


gy 


L 82 


where 


r * 

gx 


"Rcos(X)cos(GHA - 6)‘ 

r ’ 

gy 

= 

Rcos(X)sin(GHA - g) 

r* 


Rsin(X) 

- 


- 


The elements of the precession matrix, given in Table 7, are f '• 
Ref. 21. An approximation to the nutation matrix has been used, wiiv * 
6f and 6e are taken from the ephemeris DE 96, while e is computed as in 
Table 6. 

The approximations to the precession and nutation matrices are 

9 

accurate to a few parts in 10 in both angular position and radial 

magnitude. Since |r | 6 x 10 km, the corresponding error in arrival 

g 

time is on the order of picoseconds. 
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Table 7. Elements of the Precession Matrix 


Pll = 1. 

0 - T^(2.9697 X lO"'^ + 1.3 x 10“^T) 


r -2 -6 "6 1 

Pl2 = T 

2.234988 x ID + T(6.76 x 10 - 2.21 x 10 T) 

Pl3 =T 

9.71711 X lO"^ - T(2.07 x lO”^ + 9.6 x 10“^T)j 

P = -P 

21 

12 

P22 = 1 

- T^(2.4976 X 10"^ + 1.5 x lO"^ T) 

?23 = -T^ (1.085 X 10"^ + 3 X 10 ^T) 

P = ^P 

‘31 

13 

P = P 


32 ‘23 


2 -5 -8 

"33 = 1 

- T (4.721 X 10 - 2 X 10 T) 


Note: T is in Julian centuries of 36525 days past 2433282.423. 

Larger errors appear when referring the arrival times to the solar 
system barycenter. Since one normally works in 1950.0 coordinates and 
solves for pulsar positions and proper motions in these coordinates , 
barycentric arrival times are not highly sensitive to the recent increase 
in the accepted value of the general precession in longitude by 1.1 arc-sec/ 
century. However, changes in the Earth’s barycentric position can vary 
significantly as ephemerides are Improved. For example, the difference 
in the angular position of Earth in DE 108 and that in DE 96 adds up to 
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about 0.5 arc-sec/century.* In a ten-year observing program, the differ 
ence is then 0.05 arc-sec, similar to the errors in pulsar positions 
derived from the arrival time data. 


C. DISPERSION 


Considerable effort was expended in the early observations to 
obtain one frequency, 2388 MHz. This was not always possible, for 
several observations occurred at 2295 and 2328 MHz. All observations 
are now at 2295 MHz. Integration of Eq. (1) provides the corrections 
in time to be added to an arrival time measured at F MHz to refer it 
back to 2388 MHz. Performing the integration, the correction becomes 


ATp = 


DM 

2.410086 X 10 


-4 


sec 


(2388)' 


(19) 


where the values of DM are from Table 4. 

In the application of AT^ to PSR 0833-45 it was found that an 
add. '.ional 250 psec had to be subtracted frbm the measurement when 
F = 2295 MHz. No significant addendum was needed when F = 2328 MHz. 

The average pulse shape must be independent of observing frequency if 
the correction AT^ is to be perfectly adequate. Most shapes do change 
significantly with frequency, including that of PSR 0833-45, thus 
creating the need to measure an addendum to ATj^. 

D. DOPPLER DISPERSION 

The actual frequency of observation varies with the component of 
antenna velocity laying along the direction of the pulsar. The measured 
arrival time then has a varying component since the Interstellar medium 


*E.M. Standlsh, private communication. 
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is dispersive. The Doppler-shifted observing frequency is related to 
the mean observing frequency by 



where is the barycentric velocity of the antenna. In this situation, 
integration of Eq. 1 leads to a form similar to Eq. (19). Performing 
that operation. 


AT 


DD 


DM 

1.205043 X lO'"^ Fq 



1 


where fhe correction for Doppler dispersion is added to the meas- 

ured arrival time. Using the fact that |V^/c| < 10 , and setting F^ = 

2388 MHz, the correction becomes 

ATpp » 1.455220 x 10“^ ) /c (20) 

The magnitude of this yearly oscillation is about 6 psec for 
PSR 0833-45 and about 18 psec for PSR 1933+16. This correction is not 
small in all cases compared to the measurement accuracy, nor is it small 
compared to the desired uniformity of the approximation to ephemeris 
time. 

E. ARRIVAL TIME IN UTC 

The epoch at which the data sampling process begins is defined by 
a 1-sec pulse from the station clock. A sequence of 1-secoud pulses 
(clock pulses) is sent continuously via microwave link from the Standards 
Laboratory at the Goldstone complex to the individual tracking stations. 
Allowing for the delay over the links, these pulses define a UiC which 
is within 5 psec of UTC as disseminated by the National Rureat of 


68 



Standards (NBS) . The 1-second pulses generated at DSS 14, UT(DSS 14), 
are kept within 5 ysec of the Standards Laboratory pulses. Therefore, 
at DSS 14 the 1-sec pulses are assumed to be within 10 psec of UTC (NBS) . 
The 1-sec pulse used at DSS 13, UT(DSS 13), is frequently synchronized 
with the Standards Laboratory pulses to within 1 psec. A piece-wise linear 
curve approximating the deviation in this 1-sec pulse from UTC (NBS) is 
used to correct epochs at DSS 13 to UTC (NBS) . The correction is at 
most 100 psec. 

The correction equal to the width of the clock pulse at DSS 14 is 

subtracted as part of this correction. (See the discussion of system 

codes 2, 3, and 5, and Fig. 3.) The small correction equal to P /5000 

a 

is then added if appropriate. 

F. CONVERSION FROM COORDINATED UNIVERSAL TIME TO EPHEMERIS TIME 

Coordinated Universal Time (UTC) is an approximation to a uniform 
time based on the rotation of Earth. Several steps are necessary to con- 
vert to the truly uniform ephemeris time (ET) . Of prime importance 
is the correction for annual effects due to the orbital motion of the 
Earth-based clocks about the Sun. The correction AET if given in 
(Ref. 16) as 

AET = 0.001658(oin E + 0.0368) sec (21) 

where E is the eccentric anomaly of the Earth in a heliocentric orbit. 
Letting M represent the mean anomaly of Earth and e the orbital 
eccentricity, 

sin E «• sin M (1 h e cos M) 


69 



Letting an ET epoch be represented by a Julian ephemeris date (JED) plus 
the number of seconds (SEC) past 0^ ET, 

M * 358.000682 + 0.9856002628d - T^(1.55 x 10‘^ +3.33x 10“S) degrees 

-5 _7 

e = 0.0167301085 - T(4.1926 x 10 + 1.26 x 10 T) 

where 

d = (JED - 2433282.5 + SEC/86400) days 
T = d/36525 centuries 
(see Ref. 17) 

Eq. (21) omits a monthly term of 1.5 ysec, a 12-year term of 21 ysec, 
and a 29-year term of 5 ysec (Ref. 10). 

The ET epoch corresponding to each UTC arrival time, represented by 
a Julian date JD plus the number of seconds SEC past 0^, is computed 
from 


ET - UTC(USNO) = 


AET+ 38.662738 + 0.002592d sec; JD < 2441317.5 

( 22 ) 


V AET + 32.184 + n sec; JD > 2441317.5 
where AET is computed using Eq. (21), d = JD - 2440000.5 + SEC/86400, n 


is an integral number of seconds, increasing once or twice a year by 1 sec- 
ond (Ref. 15), and UTC(USNO) is the UTC kept by the U.S. Naval Observa- 
tory. The following assumptions and relations were used in computing 
Eq. (22) before 1 Jan. 1972! 

(1) lUTC(NBS) - UTC(USN0)| < 10 ysec (Refs. 12, 13) 

(2) A.l - UTC(USNO) » 6.5131177 + 0.002592d sec (Ref. 14) 

(3) lAT - A.l = -0.03439 sec (Ref. 15) 

(4) ET - lAT * 32.184 sec (Ref. 15) 
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After Jan. 1972, the appropriate number of Integral seconds Is added to 
the UTC epoch to obtain lAT (Ref. 15) under the assumption that |UTC(NBS) 
- UTC (Bureau de l'Heure)| < 10 ysec (Ref. 12). Then ET is obtained from 
ET - lAT = 32.184 sec. 

G. AN EQUIPMENT IDIOSYNCRASY 

The equipment problem of system code 3 requiring a correction of 
-0.996550 seconds is taken into account along with the corrections 
discussed above. 

H. A TEST OF THE MEASUREMENT CONSISTENCY 

The consistency of the measurements is determined by comparing the 
results using the procedures described earlier in the editing process 
(Sec. III). This comparison led to the correction of 0.996550 sec for 
the equipment failure, and the additional 250 ysec correction at 
2295 MHz for PSR 0833-45. The question finally arises concerning the 
consistency of arrival times measured at different stations. Six epochs 
exist for PSR 0833-45 for which measurements were made at DSS 13 and 
DSS 14 within 1 day of each other. Results of three of these measure- 
ments showed DSS 14 results differing signiflcat.tly from those at DSS 13. 
There were, however, serious receiver tuning problems at DSS 14 at two 
of the three epochs. Of the three results showing agreement, two are 
presented in Figure 37, where the difference (residual) .between the 
measured and predicted arrival times is presented for each Julian date 
(JD) of observation (MJD ■ JD - 2440000.5). The agreement is remarkably 
good, leading to the conclusion that no systematic differences between 
DSS 13 and DSS 14 measurements remain in this data set, 
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Figure 37. Time-of-arrival Residuals of Measurements Made at 

DSS-13 arid DSS-14 on the Modified Julian Date (MJD) 
of Observation for PSR 0833-45 




SECTION VII 


THE TABULAR RESULTS 

An example of the observational results is presented In Table 8. 
Columns 1 and 2 list the ephemeris time of the pulse arrival at the 
geocenter as a Julian ephemeris date and the number of seconds past 
0 (ET) . (Note that no complete removal of the dispersion effect has 
been applied. That is, these arrival times refer to F - 2388 MHz, not 
F = CO.) The estimate of the expected standard deviation of the measure- 
ment (defining the 67 percent confidence interval if the measurement 
errors are Gaussian) appears in Column 3. In Column 4, all of the time 
conversions from UTC to ET discussed in Section VI are listed. Column 5 
contains the applied drift correction, while the configuration code 
appears in Column 6. Barycentric arrival times are listed in Column 7 
for analysts satisfied with the accuracy of DE 96. 
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